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1 


INTROOUCTION 


1.1  Numerical  methods  for  linear  boundary  value  problems  for  CMrdlnary 

differential  equations 

The  standard  techniques  for  the  numerical  solution  of  BVP's  for 
ODE'S  can  be  divided  into  two  classes.  On  the  one  hand  are  the  "direct" 
methods  based  on  various  versions  of  finite  differences,  finite  elements, 
or  collocation,  and  on  the  other  are  the  "indirect"  or  sequential  methods 
based  on  the  numerical  solution  of  auxiliary  initial  value  problems. 
Typical  of  this  class  are  various  shooting  and  multi-shooting  approaches. 

Direct  methods  are  characterized  by  the  solution  of  global  (linear) 
algebraic  systems  for  the  discrete  solution,  (in  this  sense  multi¬ 
shooting  may  be  regarded  as  a  hybrid  method  between  the  two  classes.  A 
sophisticated  recent  example  of  multi-shooting  is  the  BOUNDPAC  package^of 
Mattheij  and  Staarink  [17])  The  solution  of  the  linear  systems,  which 
should  of  course  be  regarded  as  part  of  the  numerical  solution  of  the 
original  problem,  can  be  accomplished  by  an  iterative  method  or  by  a 
direct  method  based  on  elimination.  The  discretization  itself  may  be 
adaptive  or  non-adaptive ;  adaptivity,  however,  generally  requires  multiple 
solutions  of  the  problem.  The  well-known  programs  COLSYS  [1]  and  PASVA 
[16]  are  based  on  methods  from  this  class. 

'Indirect  solution  methods  are  characterized  by  the  association  of 
the ^BVP^  with  certain  auxiliary  initial  value  problems  (IVP).  The 
auxiliary  IVP's  are  generally  solved  uni-directionally  (forward),  but  a 
subclass  of  initial  value  based  methods  are  based  on  bi-directional 
(double  sweep)  strategies.  We  use  the  term  factorization  for  this  sub¬ 
class,  and  it  is  to  the  analysis  and  exploitation  of  the  underlying 
structure  of  the  methods  of  this  class  that  this  paper  is  addressed. 


Both  classes  mentioned  above  are  well-represented  in  the  litera¬ 
ture.  In  [2],  [3],  [13]  and  [15]  the  relationship  between  methods  of  the 
two  classes  is  explored.  The  connection  is  made  through  the  notion  of  the 
closure  of  a  numerical  algorithm  ([5],  [6]).  It  can  be  shown  that  if  one 
incorporates  the  numerical  solution  of  the  discretized  (algebraic)  equa¬ 
tions  into  the  algorithm  for  direct  methods,  then  the  algorithm  can  be 
interpreted  as  the  application  of  a  special  sequential  numerical  solver 
for  some  naturally  associated  initial  value  problems.  We  refer  here  to 
[3].  [11],  C13],  [15]  and  [18]  for  some  recent  papers  addressing  in 
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various  ways  this  relationship.  In  this  paper  we  explore  the  subclass  of 


indirect  methods  called  factorizations,  discuss  the  principles  of  their 
adaptive  construction  by  the  computer  itself,  and  address  certain 
questions  of  implementation.  Numerical  examples  illustrate  the 
effectiveness  of  a  general  code  based  on  the  method. 


1 .2  The  BVP  and  the  goals  of  the  computations 

Consider  the  linear  two  point  boundary  value  problem 


s.  S  S  ^  So, 


w’(a)  *  B(3)w(3)  -  F(s), 


Uiw(s^)  -  u^ ,  U2w(s2)  •  U2, 


where  B(*)  is  an  n  *  n  matrix  function,  w(*)  and  F(*)  are  n-vector 
functions,  and  U2  are  n.,  *  n  and  02  *  n  matrices,  respectively, 

and  u^  and  U2  are  and  n^-  and  n2-vector3,  respectively. 

The  goals  of  the  computations  are  as  follows:  given  a  set  of 
target  points 


s.  ”  a.  <  a.  <  •  •  •  <  a 


w*  V 


and  a  tolerance,  t,  find  vectors  Wj^  such  that 

A 

(3)  Wj  *  w^(o^),  i  -  0,...,n 

where  w^(*)  is  the  exact  solution  of  a  perturbed  problem 

(4a)  w^(s)  *  (B(s)  +  bj^(s))Wj^(s)  -  (F(3)  +  fi(s)),  3.|  S  s  ^  32, 

(4b)  (U  +  V  )w  (s  )  -  (u  +  V  ),  a  -  1,2, 

a  a  1  a  a  a 

where  the  perturbations  may  depend  on  the  target  point,  as  indicated,  but 
satisfy 

(5)  |b^(.)|,  lf^(-)|,  |VJ,  IvJ  ^  T. 

The  norms  in  (5)  are  a-priori  selected  and  may,  in  the  case  of  bj^  and 
f^,  be  of  the  type  Lp,  1  ^  p  S  «.  The  vector  Wj^  is  a  trace  of  an 

exact  solution  of  the  problem  (1)  with  perturbed  input  data;  the  perturba¬ 
tions  depend  on  and  may  be  different  for  different  i. 

Obviously  the  aim  of  the  computation  is  directly  related  to  the 
interpretation  of  the  numerical  solution  of  engineering  problems  where  the 
input  data  are  not  known  precisely.  That  is,  the  class  of  perturbed 
problems  (4)  are,  for  perturbations  of  a  known  magnitude  t,  indistin¬ 
guishable  with  respect  to  the  engineering  interpretaion  of  their  exact 
solutions. 

1.3  Factorization  methods  and  their  adaptive  construction 

We  consider  a  class  of  methods  based  on  the  bi-directional  adap¬ 
tive  solution  of  IVP's  for  certain  ODE’s  which  themselves  are  selected 
adaptively.  If  these  equations  were  solved  exactly  then  the  exact 


solution  of  (1)  would  be  obtained  at  each  target  point.  Let  us  assume 
that,  In  the  forward  direction,  say,  these  IVP’s  are  of  the  form 

(6a)  <t'(s)  -  R(3,<t(3)),  3^  i  s, 

(6b)  4(3^)  ■  . 

We  cannot  obtain  <t  exactly,  but  only  an  approximate  solution  $  which 
solves  (exactly)  the  equation 

(7)  5'(s)  -  R(3,i(s))  +  r(3) 

where  |r(*)|  2  t,  the  given  solution  tolerance,  and  j*!  is  a  suit¬ 
able  Lp  norm.  In  [7]  It  is  shown  how  various  local  error  control 
strategies  achieve  |r|  S  t  for  different  Lp  norms  using  a  minimal 
number  of  steps. 

We  study  a  class  of  methods  which  directly  ties  the  norm  |r|  to 
the  perturbations  of  the  input  data  of  the  original  BVP  in  the  sense  that 
bih  IfJ  Ivj.lvj  ^  C|r|  where  C  is  an  a-priori  known  constant  inde¬ 
pendent  of  the  problem  (1)  with  C  »  1.  Not  only  are  the  ODE's  (7)  solved 
adaptively  in  order  to  ensure  the  tolerance  t,  but  the  ODE’s  themselves 
are  constructed  adaptively  in  order  to  ensure  that  C  »  1 .  It  is  in  this 
precise  sense  that  the  methods  we  consider  are  optimal. 

1.4  Outline  of  the  paper 

In  Section  2  we  formulate  and  analyze  the  class  of  methods  which 
satisfy  the  requirements  stated  above.  Section  3  focuses  on  the  sub-class 
of  factorizations  which  reduce  to  matrix  Riccati  equations.  We  derive  and 
analyze  a  special  solver  for  such  matrix  Riccati  equations.  The  Riccati 


solver  is  the  foundation  of  a  general  linear  TPBVP  code  whose  performance 
is  Illustrated  on  a  varied  selection  of  example  problems  in  Section  4. 

The  efficient  and  cost-effective  solution  of  linear  problems  has 
enabled  this  method  to  be  applied  also  to  large-scale  nonlinear  BVP's  with 
turning  points  and  bifurcations.  The  description  of  this  approach,  its 
analysis,  and  our  computational  experience  will  be  reported  elsewhere. 


2 


THE  FACTORIZATION  METHOD 


2.1  The  linear  two-point  boundary  value  problem 
Definition  1 . 

(a)  Matrix  B  is  said  to  be  of  size  (m,n)  if  it  has  m  rows  and 
n  columns. 

(b)  Matrix  B  of  size  (m,n)  is  said  to  be  of  type  (m,n)  if  it  has 

maximal  rank.  9 

Suppose  that  bounded  and  measurable  matrix  functions  s  B(s)  of 
size  (m,n)  and  s  -»  F(3)  of  size  (n,1),  n  >  0,  are  given  on  the 
interval  Csi,S2]  ^  •  Suppose  also  that  for  o  -  1,2  matrices  of 

type  (n^j.n)  and  u^  of  size  with  n^  +  n2  =  n  are  given.  We 

seek  an  absolutely  continuous  (a.c.)  vector  function  s  w(s)  of  size 
(n,1)  on  [3^,32]  which  solves  the  first  order  linear  two  point  boundary 
value  problem  (TPBVP) 

w'(3)  =  B(3)w(3)  -  FCs),  a.e.  s  (  [3^,32!, 

(1) 

U^w(s^)  •  ,  U2w(32)  “  U2. 

We  will  refer  to  the  TPBVP  (1)  as  P(B,  FjO^^.u^^.s^j) . 

The  separated  boundary  conditions  in  ( 1 )  are  no  real  restriction, 
for  the  mixed  problem 

1w'(s)  -  B(s)w(3)  -  F(3),  a.e.  s  €  [3^,32! 

U^w(3^ )  +  U^w(3^)  =  u 

can  easily  be  cast  into  the  form  (1).  The  algorithms  for  the  solution 
of  P  that  we  consider  are  based  on  the  integration  of  certain  associated 
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initial  value  problems.  Accordingly,  we  turn  to  the  development  of  some 
auxiliary  results  concerning  matrix  initial  value  problems. 

2.2  Auxiliary  results  for  matrix  initial  value  problems 

Lemma  1.  Let  matrix  functions  s  ♦  3(s),F(s)  of  size  (n,n)  and  size 
(n,1),  respectively,  be  given  on  [3^,32].  Suppose  that  matrix  U  is  of 
type  (p,n),  p  S  n  and  that  u  is  of  size  (p,1).  Suppose  further  that 
s  -*■  v(s)  of  size  (n,1)  is  absolutely  continuous  and  satisfies 


■  •.■"■7 

•'■v-y-! 


*  .  •  .  < 


(1) 


v'(3)  =  B(s)v(3)  -  FCs) 


a.e.  3  €  ,  r^]  [s^  ,3^] 


UvCjq)  =  u  for  some  €  [0^,02]. 


If  3  -►  $(3)  of  size  (p,n)  defined  on  lOyCi^l  is  a.o.  and  satisfies 


(2) 


|$<(3)  =  -4'(s)B(3)  +  Z(3)i>(s),  a.e.  s  €  la^,o^2 


^^COq)  =  U, 


?-v 


.4 
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y:-v 

ii  .%  .N . 


and  if  3  ^  cpCs)  of  size  (p,1)  is  a.e.  and  satisfies 


(p'Cs)  =■  -<J(3)F(3)  +  Z(s)<P(s),  a  e.  s  €  [0^,0^] 


(3) 


where  s  Z(s)  of  size  (p.p)  on  [0^,02^  continuous,  but  otherwise 
arbitrary,  then. 


(4) 


<K3)v(3)  =»  (P(3) 


Vs  €  [(j^  ,  02] 


Proof.  Define  \p(s)  on  [3^,02]  by 


Then 


and 


41(3)  »  (<frv  -(p  )(s) . 

ii;(  Jq)  =  Uv(  <Jq)  -  u  *  0, 

.  $'v  +  'jv’  —  (pf 

=  (-$8  +  Z<t)v  +  $(Bv  -  F)  -  ( iPF  +  Z  ) 


-  Z(<tiv  -  (p)  »  Zip,  a.e.  3  € 

Since  Z  is  continuous,  ip  *  0  by  uniqueness  and  the  lemma  is  proved. 

The  method  of  factorization  for  BVP’s  is  based  on  the  propagation 
of  the  boundary  condition  across  the  interval  in  a  manner  consistent  with 
the  differential  equation  (1);  Lemma  1  is  the  formal  expression  of  the 
nature  of  this  propagation.  We  use  the  following  terminology . 

Definition  1.  A  matrix  function  $  of  size  (p,n)  satisfying  (2)  is  a 
transition  matrix  based  at  ^Oq.  The  matrix  Z  which  induces  ip  is  the 
associated  conditioning  matrix  on  [0^,0^]  ^  C3^,S2]-  The  vector  func¬ 
tion  cp  of  size  (p,1)  satisfying  (3)  is  a  transition  vector  based  at 

Cq.  U 


Lemma  2. 

Suppose  that  <t 

and  f 

are  size  (p,n)  transition  matrices  on 

based  at  Oq,  Z 

and  Y 

are  continuous  size  (p,p)  condition- 

ing  matrices  on 

that 

(5) 

<p  *  s 

-<PB  Zip, 

(6) 

41  '  s 

-fB  +  Y'y, 

8 


on  [a^jO^],  and  that 

(7)  -  K„»(o„) 


with  Kq  of  type  (p,p).  Then  there  exists  a  matrix  function  s  -*•  K(s) 
of  type  (p,p)  on  [0^,02]  such  that. 


(8) 


$(3)  =  K(3)  'P(s), 


s  €  [  0^  ,  02  ] 


Proof.  Let  s  KCs)  solve  the  (linear)  initial  value  problem 


(9) 


K'(3  )  =  Z(3)K(s)  -  K(s)Y(s), 


K(ao)  =  Kq. 


Then  we  have  that  on  [0^,02] 


(10)  (KH'-*)  ’  -  -(K'i'-'t>)B  +  Z(K'i'-'t) 

with  initial  condition 


(11)  (:<¥-<t)(0Q)  =  0. 

By  uniqueness  of  solutions  to  (10),  (11)  we  have  that 

(12)  Hs)  =•  K(3)'I'(3),  3  €  [0^,02] 

It  remains  to  show  that  K  is  of  type  (p,p),  i.e.,  invertible.  Let 

L(s)  be  the  solution  of 


L'(3) 


Y(3)L(s)  -  L(3)Z(s)  , 


then  KL  satisfies 


(13)  (KL)'  -  Z(KL)  -  (KL)Z 

(14)  (KL)(<Tq)  -  Idp. 

Since  the  p  *  p  identity  matrix  Idp  is  the  unique  solution  of  (13), 

(14)  on  [0^,02],  the  result  is  proved.  0 

Leona  3.  The  rank  of  *(s)  satisfying  (5)  is  constant  on  [0^,02]. 
Proof.  Let  3  -*■  E(3)  of  size  (n,n)  be  the  solution  of 

E'(3)  =  -E(3)B(s)  3  €  [0^,02] 

E(ao)  -  Id^. 

Clearly  't  solving  (6)  with  Y  =  0  exists  and  is  given  by 

^(s)  -  '1'((7q)E(s), 

so  that  by  Lemma  2 

(15)  “Ka)  =  K(3)H'(0q)E(3). 

If  we  show  that  E  is  non-singular  for  s  €  [0^,02]  then  from  (15)  we 
can  see  that 

rank  9(3)  -  rank  ^(s) 

which  proves  the  result. 

To  show  that  E(3)  is  of  type  (n,n),  let  F(3)  solve 


(16) 


F'(3)  -  B(s)F(s) 

F(ao)  -  Id^. 

10 
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It  is  easy  to  see  that  E(3)F(3)  -  on  [0^,02]. 


We  apply  Lemma  3  to  show  that  solvability  of  (4)  at  one  point 
Oq  €  [0^,02]  gives  solvability  at  each  s  C  [0^,02].  In  fact  we  can 
prove 

Corollary  4.  For  <t,  <j>  satisfying  (2),  (3),  if  the  equation 


UWq  -  u 

has  k  independent  solutions,  then 

♦(s)w  »  (p(s) 


has  k  independent  solutions  for  every  s  € 

Proof.  Apply  Lemmas  1  and  3  to  the  augmented  matrix  which  solves 


c*:cp] 


(17) 


-C  (p] 

I 


F 

OJ 


[  «  1  <p]  (  df 


culu]. 


By  by  hypothesis,  rank  U  =  n  -  k  -  p  and  u  €  range  U.  Now 

rank[<P  (P]  -  p  and  rank  $  =  p  by  Lemma  3.  Thus  tp  €  range  $  and  the 

proof  is  complete. 


2.3  Definition  of  a  factorization 

Let  boundary  value  problem  P(B, F,U(j,Ujj,Sg)  be  given  as  in  Section 

2.1 


%  *. 

•  •  fc 

V,* 

p  V  V 


Definition  1.  A  factorization  of  P  consists  of; 


partitions  ir  ,  a  -  1,2  of  [3^,82]  with 


ir^  :  3^ 


(0)  ,  (1)  ,  ,  ^“1^ 

0  v***x(J  *3j 


'a- 


'"2>  ,  'V>  ,  ,  (0) 

O2  ^  *^2  <  •  •  •  <  (J2 


^2’ 


collections  Z^,  a  »  1,2  of  conditioning  matrices  of  size 


(n  ,n  ), 
a  o 


{3  -  Z^^^3):  3  €  i  -  1  ,m^}. 


where  we  use  the  notation 


r  (i-1)  (i),  . 

t  J,  a  ■  1  , 


r(i) 


.  (i)  „(i-1)T  „  ,  .. 

L  cfg  » (^2  >  cx  ■  2  , 


collections  ^  ,  a  »  1,2  of  transition  matrices  of  type  (n  ,n) 
01  a 


3  >  i  -  1,m  }; 

a  a  a  a 


collections  <P  ,  a  -  1,2  of  transition  vectors  of  size 

a 


■  (s  -•  q>“’(3);  3  «  I*‘*,  1  - 

a  cx  a  a 


,  a  =  1 ,2  of  scaling  matrices  of  type  (n  ,n  ) , 


collections  K 


(Ff)  collections  P  ,  a  -  1,2  of  constant  similarity  matrices  of 

Cl  ‘  ^  '  ■  ' ' 


type  (n,n), 


1  -  I,.,), 


such  that  the  following  conditions  hold  for  i  =  Itm^^  and  a  » 


B^^^Cs)  -  (P^^^)“^B(s)P^^\ 

Cl  a  a 


s  €  li 


a 


(P^^^“^F(3),  3 

a  ® 


'(s) 

a 


-<Ii^^^(3)B^^\3)  +  Z^^^(s)il>^^\3), 
a  a  a  a 


3  € 


a  a  a  a 


3  *  I<‘> 


a  a  a  a  a  a 


(i),  (i-l),  „(i-1), (i'1),  (i-1). 

iFo)  (D  (0  j-K  <p  (a 
^a  a  a  a  a 


In  order  for  (F5)  and  (F6)  to  make  sense  for  i  =  1,  we  have  used  the 
notational  conventions 


(F7) 

a  a 


(0),  (0).  , 

(p  (a  )  »  u  . 
a  a  a 


The  initial  value  problems  (F3),  (F5)  and  (F4),  (F6)  for  $' 


and  respectively,  are  posed  forward  in  s  for  a  ■  1  and 

backward  in  s  for  a  -  2.  We  combine  the  forward  and  backward 
factorizations  into  a  composite  factorization  as  follows; 

Let  the  composite  partition  ir  -  ir^  U  1T2  be  given  by 


(1) 


.(0) 


.(1) 


,(m) 


-  s 


2’ 


and  let  i  *  1  ,m.  It  is  obvious  that 

Lenna  1.  There  exist  unique  indices  i^,  a  «  1,2  such  that 


r(i) 


n  I 


(i,) 


For  a  €  ^  ^2  define  '(0)  of  size  (n,n)  by 


(2) 


.(i) 


(i.)  (i,)  -1 

4,  (a)(P^  ) 


(i-)  (ip)  -1 

«2  (o)(P2  ) 


<p^^^(o)  of  size  (n,1 )  by 


and  then  set 


(6c) 


Finally,  we  set 

(7) 


C 

o 

0 

i 

n^,n^ 

1 

1 

1 

and  then  state 

Lemma  2.  is  of  size  (n,n),  is  of  size  (n,1),  is  oi 


size 

(n,n),  is 

of  type  (n,n),  R^^^  is 

of  type 

(n,n) 

(and  so 

K(i) 

is  well-defined). 

and  we  have 

(8) 

- 

V  ■ 

Al)a  *  7(i)*(i)  o  a 
-  B  +  Z  V  a.e. 

on 

i  - 

(9) 

ep  « 

-  ■P^^^F  +  a.e. 

on 

i  - 

1 

(10) 

^i-1^(i-1)(^(i-1))^ 

i  • 

1  •  %  •  •  1  ni^ 

(11) 

,l-1^(l-1)(,(i-1)). 

i  - 

1  f  •  •  •  y  tQ^ 

Proof 

.  Equations  (8)-(11)  follow  directly  from 

(FI  )-(F6). 

We  also  establish  a  notation  for  the  composite  factorization. 
Definition  2.  A  composite  factorization,  F  for  P(B, F,Ujj,Ujj^,3|j)  con¬ 


sists  of  the  sets 


,  {g  +  $'■^^(3)  :  3  satisfies  (8)}, 

(P  -  [3  ♦  (P^^^s)  :  3  6  <p^^^  satisfies  (10)}, 

Z  -  {3  >  Z^^^Cs)  :  3  Z^^^  as  in  (4)}, 

K  =>  as  in  (7)}. 

We  write  F-  F(  ir,  *,  <P,  Z,K)  *  {ir,<t,4>,  Z,K}.  p 

Theorem  1-  Suppose  F(ir,iP,{p,Z,K)  is  a  factorization  for 
p( B, F,U^,u^,3jj)  and  let  3  ■*■  wCs)  be  an  a.c.  solution  of  P.  Then 

(12)  $w  »  {p. 

The  sense  in  which  (11)  is  to  be  interpreted  will  be  apparent  from  the 

Proof.  Let  a  target  point  a  €  [s^,S2]  be  given.  Then  a  f  foi 

some  index  i,  0  S  i  S  m.  We  must  show  that 


(13) 


<P^^^a)w(o)  -  (P^^^a). 


Since 


(i^) 

n  it  is  enough  to  show  that 


(i^)  (i„)  -1  (i  )  (i  ) 

(14)  <S  “  (3)(P  “  )  w(3)  -  cp  “  (3)  Vs  €  I  “  ,  a  =  1,2. 

a  a  a  a 


We  prove  (14)  by  induction  on  i^^. 

i(j  -  1  :  By  definition  of  P  we  have  that 
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and  so 


U  w(o^°^) 
at  at 


K«>'UW(,<°>) 
a  01  a 


Cl  a 


But  by  (F5)  and  (F7)  we  have 


Cl 


^(0)  p(1) 

01  01  a 


and  by  (F6)  and  (F8) 


CL  CL 


a  a 


Now  by  virtue  of  Lemma  2.2-1,  (111)  holds  for  =  1 .  In  particular,  (14) 
holds  for  s  =  and  la  ’  1  •  Since  ^  is  of  type  ("a’^o^’ 

gives 

(15) 

a  a  a  a  a  ^  a 


^01  ^  *  Having 


(i  )  (i  )  (i  )  -1  (i„)  (i„)  (i^) 

K  “  ♦  ®  (P  “  )  w(o  “  )  -  K  “  4>(a  “  ) 

CL  CL  CL  CL  CL  CL 


and  (F3)-(f6)  for  i  -  i^  +  1  we  again  apply  Lemma  2.2-1  to  conclude  that 
(14)  holds  for  i  =  i(j  +  1 .  b 


The  proof  of  Theorem  1  actually  gives  us  a  bit  more  than  is  stated 

in  the  theorem.  The  induction  together  with  the  observation  that  each 

K  is  of  type  (n  ,n  ),  i.e.,  invertible,  and  an  appeal  to  Lemma 

a  a  a 

2.2-3  gives 
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Cwollary  2.  The  rank  of  a  »  1,2,  is  constant  on  [3^,32].  That 

is,  'Sjj  i3  of  type  {n^,n)  on  [31,32]. 

Now  that  we  have  shown  that  (11)  holds  if  w  solves  P,  we  show 
that  the  solvability  of  (11)  is  also  inherited  from  P.  In  fact,  they  are 
equivalent. 

Theorem  3-  Let  F(7r,'t,<p,Z,K)  be  a  factorization  for  P(B,  F,U^,u^,s^) . 
Then  for  0  €  [si,S2]i  the  linear  algebraic  equation 

(15)  <|i(a)w  =  <p(o) 

has  as  many  linearly  independent  solutions  w  €  R"  as  there  are  linearly 
independent  a.c.  functions  s  -*■  w(3)  solving  P. 

Proof.  Equation  (15)  is  interpreted  in  the  sense  of  Theorem  1. 

If  s  w(s)  solves  P  then  by  Theorem  1  w  =  w(o)  is  a  solution 
of  (15). 

On  the  other  hand,  let  Wi  be  a  solution  of 

'S(s^)w,  »  <J>(s^) 

and  let  s  -*■  w(s)  be  the  solution  of  the  initial  value  problem 

Iw'(s)  “  B(s)w(s)  -  F(s),s  >  Si 

w(si)  -  Wi . 

By  the  definition  of  ^i  and  cpi  and  by  the  fact  that  is  of  type 

^■^1  »ni )  we  have 

Uiw(3i)  =  Ui . 
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We  must  show  that 


(16)  U2w(s2)  •  U2. 

Now  by  Lemma  2.2-1  applied  with  p  -  n  »  +  n2,  U  »  <t(s^),  and 

u  -(p(s^),  we  conclude  that 


from  which 


<t(s2)w(32)  •  <p(s2) 


<l>2(s2)w(s2)  -  (P2(S2)- 


Therefore 

K^'^^U2w(s2)  =  4°^U2 


and  since  us  of  type  (n2.n2)  we  have  (16).  B 

In  our  development  above  we  have  assumed  that  a  factorization 

for  P  exists.  In  fact,  there  are  many,  as  we  will  show  by  example 

below.  The  trivial  factorization  (Z  =  0,  K  =  always  exists,  for 

example;  it  is  usually  not  numerically  realizable,  for  it  is  the  shooting 

method.  A  factorization  algorithm,  then,  should  select  the  conditioning 

matrices  the  scaling  matrices  the  similarity  matrices 

a  a  a 

and  the  partitions  adaptively  in  order  to  ensure  the  numerical 

stability  of  the  computations. 


2.4  Stable  factorizations 

From  the  point  of  view  of  practical  computations,  the  boundary 

value  problem  P(3,F,U  ,u  ,s  )  is  but  a  model  of  physical  reality.  The 

a  a  OL 

data — the  arguments  of  P — are  by  definition  known  inaccurately,  and  it  is 


only  after  the  influences  of  these  inaccuracies  on  the  solution,  w, 
to  p  are  acknowledged  that  a  reasonable  interpretation  of  the  meaning  of 
even  the  exact  solution  to  P  can  be  made.  This  interpretation  is  based 
on  our  confidence  both  in  the  model  itself  and  on  the  reliability  of  the 
data  supplied.  We  take  the  point  of  view  that  the  exact  solutions  w  and 
w  of  P  and  a  perturbed  problem  P,  respectively,  are  equivalent 
if  P  and  P  are  indistinguishable  with  respect  to  the  goals  of  the 
computation. 

On  the  other  hand,  we  do  not  have  at  our  disposal  the  exact  solu¬ 
tion  w  to  P,  but  only  an  approximate  solution,  w.  However,  if  it  is 
possible  to  interpret  w  as  the  exact  solution  of  a  perturbed  problem  P, 
and  if  P  and  P  are  acceptably  close  (with  respect  to  the  goals  of  the 

A 

computation),  then  w  is  an  acceptable  approximation  to  w. 

In  order  to  formalize  this  idea,  let  us  suppose  that  a  norm  1*1^ 

is  given  on  R^.  Then  I* I  induces  a  natural  norm  I* I,  >  on 
°  '  'n  '  '(n,n) 

matrices  of  size  (n,n): 


(1) 


|A| 


(n,n) 


sup 

Ivl  =1 

'  'n 


|Av| 


We  will  need  to  measure  the  size  of  various  sub-matrices  of  A  in  a 
consistent  way.  Suppose  that  A  is  a  matrix  of  size  (p,q), 

1  ^  p,q  S  n.  Then 


(2) 


l^l(p,q)  “ 


(n,n) 


where  A  is  obtained  from  A  by  augmentation  by  zero: 


0  p 


OJ  n-p 
n-q 


course,  |y|^  .  for  v  ( 

Now  we  are  in  a  position  to  state 


Definition  1.  An  approximate  solution  w  of  P(B,F,U  u  ,s  )  is 

a  a,  a 

A 

6-acceptable  at  a  €  [si,S2]  if  w(a)  is  the  exact  solution  of  £ 

perturbed  problem  P  =  P(B,F,U  ,u  ,s  )  where 

a.  at  a 


and  the  norms  (^a)  and  (4b)  are  understood  as  norms  in  Lp,  1  S  p  S  ®. 

We  use  the  notation  6  =  {5n,5_,6,,  ,5  }. 

B  F  U  u 
a  a 

Based  on  the  notion  of  6-acceptability  of  the  approximate  solu¬ 
tion,  we  seek  to  formulate  factorization  algorithms  for  which  we  can 
relate  the  errors  introduced  by  the  algorithms  to  perturbations,  of  an  £ 
priori  specified  magnitude,  to  the  data  of  P.  The  specific  conditions 
under  which  the  perturbations  in  the  input  data  are  acceptable  depend  on 
the  problem  and  the  aims  of  the  computations,  i.e.,  on  the  choice  of 

In  fact,  our  theory  gives  more:  a  complete  characterization  of  the 
structure  of  the  perturbations  to  the  original  problem  (cf.  the  proof  of 


'  L*-  “-V, 


,*  ••  % ' 
■  y 


(4a) 

IB(-) 

A 

/•A 

(4b) 

|f(  •) 

^  S' 

.  •  •  -  • 

“  .•  *.* ' 

(4c) 

K- 

"al(n  ,1)  * 

a 

Sy  ,  a  =  1,2, 

a 

T-'TI 

(4d) 

K- 

'^al(n  ,n)  ^ 

a 

.  a  =  1,2, 
a 

.•,V. 
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Theorem  1  below).  Indeed,  the  stability  of  the  numerical  algorithm  will 
be  related  directly  to  the  stability  of  the  given  problem  P  to  perturba¬ 
tions  of  its  data. 


Definition  2.  A  factorization  ,Z,K)  for  problem  P(B,F,U  ,u  ) 

(j  a 

is  bounded  above  if  there  exist  ,  M2  >  0  such  that 


(5) 


<"•%> 


^  M  ,  a 
a 


1,2. 


T 

F  is  bounded  below  if  Is  invertible  for  every  s  f  [3^,82]  and 

there  exist  m^ ,m2  >  0  such  that 


(6) 


|(4  (a)*'^(s))“^  |.  ,  S  - 

’  a  a  '(n  ,n  )  m 

a  ct  0 


1,2, 


F  is  bounded  if  it  is  both  bounded  above  and  bounded  below.  H 

In  order  for  a  factorization  to  be  useful  for  practical  computa¬ 
tions  and  m^^  should  be  numbers  of  "reasonable"  magnitude.  If 

2  T 

M^/m^  is  large  the  condition  number  of  is  large;  effectively 

loses  rank.  and  m^  are  also  related  to  the  interpretation  of  the 

errors  associated  with  the  numerical  realization  of  the  factorization. 

Indeed,  suppose  that  a  €  [0^,0^]  is  a  target  point  in  an  interval 

over  which  we  have  the  forward  and  backward  factorizations  (a  =  1,2) 

(7a)  <l.'(s)  -  -(t  B)(3)  +  (Z  $  )(3), 

Cl  ct  ct  a 


(7b) 


«  (ff  ) 

Ct  a 


and 
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(8a) 


.  -.ii  -  i 


(P'(S)  -  -(*  F)(3)  +  (Z  cp  )(3), 

C(  a  a  Qt 


(8b) 


(Pn,^3  ) 

'*'(x  a 


Any  numerical  realizations  4'^,  of  (7)  and  (8)  may  be  viewed  as  the 
exact  solutions  of  a  perturbed  problem 


(9a) 

r(3)  > 

-Cy  B)(s)  (Z  4-  )(3)  +  A  (s), 

ct  a  a  a 

(9b) 

V  (o  )  => 

ot  a 

U  +  V  ; 
a  a 

and 

(10a) 

4»'  (a)  - 

(X 

-(4'  F)(3)  +  (Z  t|;  )(3)  +  6  (s), 
a  a  a  ct 

(10b) 

u  +  V  . 

Cl  a 

In  (9)  and  (10)  the  matrix  of  size  (n^^.n)  and  vector  6^^  of  size 
(n^^.l)  represent  the  discretization  errors  of  the  numerical  method  used 
to  solve  (7)  and  (8).  The  matrices  and  v^  represent  the  error  in 
realizing  the  boundary  conditions.  We  claim  that  s  ■*  v(s)  defined  on 
[o^.a^]  by 


’4'^(3)' 

(  3  ) 

v(3)  - 

4',(3) 

u  2 

_ip2(3). 

is  a  S-acceptable  solution  to 

(12a)  w'(3)  -  3(3)w(3)  -  F(3),  s  € 


(12b) 


U,  W(  <3,  ) 


u,  , 


U,w((Jo)  »  u~, 


with  a  5  determined  solely  by  the  conditioning  of  '\l^  and  the 


discretization  errors. 


Theorem  1.  Suppose  a  =  1,2  solving  (6)  and  (7)  satisfy 


(13a) 


,n)  ‘ 


(13b) 


|(y  (.)yT(  )-1|  <  1 

'  a  a  ’(n  ,n  )  m 


Then  v  defined  by  (11)  is  5-acceptable  (at  o  €  [0^,02^)  with 


(lHa) 


6q  S  max  {—  I A  I  /  v } , 
B  -  'm  *  a*(n  ,n) ^ 

a«i  ,d  a  ct 


(14b) 


'f  ‘  l;r  l^alcn  ,ol' 

a«t,2  a  a 


(14c) 


(IW) 


6  S  Ic.,  1^* 
u  '  a ' (n  ,  1 ) 


In  (14a)  (I4b)  the  norms  as  functions  of  s  are  assumed  to  be  of  L 


Proof.  Introduce  matrix  functions  bj^  by 


b  (3)  -  4'^(3)(’f  (3)'1''^(3))  ^A  (3); 

a  a  a  a  ct 


and  vector  functions  f„  by 


'.V 


•  \  a(Jl  mJ". 


'*  -Nk  '->  V-  > 


(16) 

f 

a 

(s) 

-  ¥j(3)(y  (s)4’J(3))“^5  (s), 

a  a  a  a 

It  is  easy  to 

see  that 

"a* 

i(/^  satisfy,  for  a  *  1,2, 

(17a) 

a 

-(¥  (B+b  ))(3)  +  (Z  4'  )(3), 

01  a  Cl  a 

(17b) 

= 

U  +  V  , 
a  a 

(18a) 

♦(3) 

= 

-('l'^(F  f  ))(s)  +  (Z  'f  )(3) 
a  01  Cl  a 

(l8b) 

w 

«  1 

u  , 
a 

on  [a^.a^]. 

Now  for 

fixed  a 

6  [a 

^,02]  let 

1 

b|(s),  3€ 

b(3;o) 

82^3)3  3  €  Co,02J» 

and 

1  f^  (s) ,  3  €  [0^ ,  a] 

f (3;o) 

1  f2(s),  3  €  [a, 02]' 

Then  V  »  v(o)  satisfying  (11)  at  3=0  is  the  value  at  a  of  the 

solution  of  the  perturbed  problem  P(B+b,  F+f,  U  +V  ,  u  +v  ,  0  ).  The 

a  a  a  a  a 

estimates  (14)  follow  easily  from  (13)t  (15)  and  (16).  | 
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Remarks . 

1.  The  perturbed  problem  P  through  which  we  interpret  the 
approximate  solution  via)  depends  on  the  target  point  o.  The  pertur¬ 
bation  is  different  at  every  target  point. 

2.  The  magnitude  of  the  perturbation  6  »  ,S  }  is 

D  r  U  U 
a  a 

independent  of  the  coefficients  of  the  original  problem,  as  long  as  a 
bound  of  the  type  (13)  holds.  In  practice,  it  is  the  adaptive  construc¬ 
tion  of  the  conditioning  matrices  which  will  guarantee  (13). 

3-  We  assume  that  the  effect  of  roundoff  errors  is  negligible.  It 
could  be  incorporated  into  the  matrices  and  and  into  the 

vectors  and  v^.  The  solution  of  the  local  systems  (11)  at  each 

target  point  also  Introduces  roundoff  errors.  We  assume  that  the  effects 
of  these  roundoff  errors  can  be  neglected  with  respect  to  the  discretiza¬ 
tion  errors  which  have  already  been  made,  especially  if  the  computations 
are  carried  out  in  double  precision.  It  is,  however,  possible  to  inter¬ 
pret  also  these  roundoff  errors  as  perturbations  of  the  input  data  of  the 
problem. 

4.  If  we  regard  the  errors  in  the  realization  of  the  boundary  data 

as  roundoff  errors  we  may  assume  that  =«  0,  ”  0  in  (9)  and  (10). 

5.  The  expressions  for  b^  (15)  and  f^  (16)  give  the  complete 
structure  of  the  perturbations.  This  is  in  fact  a  stronger  result  than 
the  claim  of  Theorem  1 . 

This,  then,  is  what  we  mean  by  a  stable  factorization:  a  bounded 
factorization  for  which  m^  and  are  of  such  a^  priori  known 

magnitudes  that  the  approximate  solution  is  5-acceptable  when  the  error 
tolerances  of  the  numerical  integrator  for  the  factorization  matrix/vector 


initial  value  problems  are  roughly  on  the  order  of  5. 


2.5  Examples  of  stable  factorizations 


2.5.1  Continuous  orthonormalization.  Since  In  view  of  Theorem  2.4-1  we 

T 

seek  to  control  the  behavior  of  we  study  the  properties  of  this 

T 

product.  Let  Q^Cs)  «  <t'(s)4  (s)  and  observe  that 
a  a  a 

(1)  Q„(3)  +  qJ(s)  -  ($  *’^)'(3). 

a  a  a  a 

Now  since 

(2)  Q  (a)  =  -<t  +  Z 

a  a  a  a  a  a 


and  4  is  of  type 
a 

(3) 


(n^j^.n),  we  can  solve  for  Z^! 


T  T  -1 

Z  -  ('^B*  +  Q)($'t) 

a  a  a  a  a  a 


Then  it  is  easy  to  prove 


1.  Let  Z^(s)  be  determined  by  (3)  for  s  -*•  Q^Cs)  bounded, 
measurable,  and  anti-symmetric.  Then  s  ■»  satisfying 

(4a)  I*'  (s)  -  (s)B(3)  +  Z  (s)*  (s) 

d  Ql  Ct  Cl 


(4b) 


t  (s  ) 
a  a 


U 

a 


has  the  property  that 

(5) 


♦  (3)<6^(3) 
a  a 


T 

U  U  . 
(X  a 


T 

We  can  assume,  without  loss  of  generality,  that  > 

T  -  y 

Indeed,  we  need  only  multiply  U  by  K  -  (U  U  )  *  .  Thus  we  see  that 

a  a  aa 

the  factorization  induced  by  the  conditioning  (3)  with  anti-symmetric 


maintains  the  constraint  that  the  rows  of  the  transition  matrix  «  are 

a 

orthonormal . 

However,  in  view  of  Theorem  2.M-1,  we  can  expect  to  realize  this 
constraint  only  approximately.  Indeed,  given  by  (3)  for  B  will  not 

do  for  the  perturbed  matrix  B  +  b^.  In  practice,  therefore,  the  approxi¬ 
mate  solution  will  drift  away  from  the  constraint  manifold,  and  provision 
must  be  made  for  periodic  discrete  re-orthonormalizations.  We  refer  to 
[19]  for  details  and  some  numerical  examples. 


2.5.2  Stabilized  continuous  orthonormalization.  This  factorization  is 

also  based  on  (1),  but  with  chosen  to  ensure  asymptotic  stability  of 

the  manifolds  *  a  U  U^. 

a  a  a  a 


Lemma  2.  Let  ^^(s)  be  determined  by  (3)  for  given  by 


(6) 


T  T 

where  <.  >  0  and  <0  <  0.  Then  the  manifold  ( ♦  ♦  >(3)  •  U  U  is 
'  a  o  a  a 

asymptotically  stable  (forward  in  s  for  a  -  1  and  backward  in  s 
for  o  -  2) . 


Proof. 


It  is  enough  to  show  that  the  solution 


f  (3)  »  U  of 
Cl  a  a 


f'(s) 

a 


2ic  (U 
a  a  a 


-  v=>> 


is  asymptotically  stable  forward  in  3  for  a  -  1 .  But  this  is  immediate 
since  the  coefficient  matrix  of  this  constant  coefficient  nonhoraogeneous 
linear  problem  is  ^ 
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Refflarks 


1.  In  an  implementation  of  the  continuous  orthonormalization  it 

T  -1  T  -1 

might  be  tempting  to  replace  (♦  ♦  )  in  (3)  by  (U  U  )  in  order  to 

oa  aa 

T 

avoid  the  expense  of  inverting  every  function  evaluation.  This 

of  course  invites  disaster,  as  has  been  noted  by  Meyer  [20  ]  in  a  similar 
context.  For  a  discussion  of  this  method  see  also  [10]. 

2.  On  the  other  hand,  it  is  possible  to  show  that  for  |ic^| 
sufficiently  large  in  (6),  an  implementation  of  stabilized  continuous 
orthonormalization  can  afford  to  commit  the  above-mentioned  crime. 
Unfortunately  having  |k^|  large  exacerbates  any  stiffness  inherent  in 
the  problem. 

2.5.3  Riccati  factorization.  The  Ricatti  factorization  is  based  upon 

partitioning  ♦  in  order  to  determine  a  maximally  nonsingular  n  x  n 
a  o  a 

submatrix.  We  will  need  some  notation  and  a  few  preliminary  results. 

For  0  <  p  <  n  we  define  matrices  Lp  and  Rp  of  type  (n,p) 
and  (n,n-p),  respectively,  by 


(7) 


It  is  easy  to  verify 


3. 


(3a) 


(3b) 


L  +  R  r”^  -  Id  , 

P  P  P  P  n 


T 

L  L 
P  P 


For  a  matrix  A  of  size  (n,m)  and  0  <  p  <  n  we  set 


Al.  -  A2^  -  rJa 

and  for  A  of  size  (m,n)  we  set 

A^  -  ALp,  A2  =  ARp. 


Evidently  if  A 

is  of  size 

(n,n) 

we  will  write 

form  as 

*1,1  *1,2 

A  ■* 

*2,1  *2,2 

where 

*1,1 

■ 

’  *1,2  “ 

*2,1 

- 

,  A2^2  ’ 

Finally,  if  0  =  ( . 0^)  is  a  permutation  of  the  integers  (1,...,n) 

we  denote  by  P®  the  corresponding  n  *  n  column  permutation  matrix; 
i.e.,  is  the  n  x  n  identity  matrix  with  columns  permuted  by  a. 

Then  if  U  is  of  size  (p,n)  we  have  that  UP®  is  a  matrix  the  columns 
of  which  are  those  of  U  permuted  by  0.  If  U  is  of  size  (n,m) 
then  (P®)^U  is  a  matrix  the  rows  of  which  are  those  of  U  permuted  by 
0.  Clearly,  (p'^)'^p‘^  -  p'^(p‘^)’^  -  Id  . 


im 


Lamia  4.  If  U  is  of  type  (p,n)  with  0  <  p  i  n,  then  there  exists  a 
permutation  o  such  that  the  entries  of  the  matrix  [(UP^)Lp]  are 

all  in  absolute  value  less  than  or  equal  to  one. 


Proof.  Since  U  is  of  type  (p,n)  there  exists  a  permutation  t  such 


*•  *• 


that  UP  L  is  invertible.  Let  a  be  such  that 
P 


IdetCUP'^Lp)]  ^  ldet(UP‘'Lp)| 


for  all  permutations  t.  Then  by  Cramer's  rule  the  elements  of 

cup'll  )  ^UP*^  have  the  form 
P 


det cup'll  ) 


detCUP  L  ) 
P 


for  various  choices  of  permutation  t. 


We  call  UP  L.p  with  a  as  in  Lemma  4  a  maximally  nonsingular  (p,p) 
submatrix  of  U. 

The  following  remarks  apply  both  to  the  a  =  1  (left  to  right) 
and  a  =  2  (right  to  left)  factorizations.  Accordingly,  we  suppress  the 
subscript  a  and  discuss  only  the  left  to  right  factorization. 

We  construct  a  solution  to  a  problem  of  the  following  type:  Given 
absolutely  continuous  functions  s  B(s),  F(3)  on  [s.,,S2]  with  B  of 

size  (n,n)  and  F  of  size  (n,l),  and  given  a  matrix  U  of  type 
(p,n)  and  vector  u  of  size  (p,1),  determine  (piecewise  absolutely 
continuous)  functions  s  -*■  $(3),  (p(3)  on  [3^,32]  such  that  4  is  of 

type  (p,n),  (p  is  of  size  (p,1)  and  such  that 


♦(s)w(s) 


((>(s) , 


3  €  Cs. .Sp] , 


for  any  function  s  -*•  wCs)  satisfying 


w'(s)  “  3(3)w(s)  -  F(3) 


[Uw(3^)  =  u. 


step  0;  Let  s^°^  =  s, ,  =■  3,  -  F,  =  U  and  let  a^''^ 

be  such  that  ^  ^Lp  is  maximally  non-singular,  can  be 

obtained  in  practice  by  performing  Gauss  elimination  with  column  pivoting 
on  U.  The  result  is  the  matrix 


P 


[Id  I 

p  I  2 


.<'>  i  »'"1 


If  the  elimination  is  actually  applied  to  the  augmented  matrix  [U  lu], 
allowing  only  pivoting  with  respect  to  the  first  n  columns,  then  we  also 
obtain,  with  v^*^^  =  u. 


)■'  v"” 
P 


Step  i  ( i  >  0 ) : 


Suppose  that  of  type  (p,n),  and  of 

size  (p,1)  are  given.  Let  be  such  that  \  )  is 


maximally  non-singular.  Compute 


^P*^  L  )  ^P° 

P 


j  ] 


[idp 


••.y.v 


Now  set  -  (P“^  )  -  (P°  )  F^^"''^ 

let  solve  the  initial  value  problems 


(15a)  $^^^'(s)  =  +  (Z^^^$^^^(3), 


(15b) 


^(l)^g(i-l))  _  y(i). 


(16a) 


<p^^^'(s)  =  -(4^^^F^^^)(3)  *  (Z^^^cp^^h(3), 


(16b) 


(i).  (i-1).  ,(i) 

cp  (s  )  »  U  , 


where  Z'^'^  is  given  by 


Z^^^s)  -  (${^^bJ^J  +  $^^^B2^J)(3), 


Recall  that  and  .  It  is  not  hard  to  see 

1  p  2  p 

that  (13)  and  (15)  imply  that 


<t|^\s)  »  Idp 


as  long  as  (15a)  holds.  Moreover,  satisfies 


(19a)  ♦^^^’(s) 


1,2  1,K2  ^2  2,2  ^2  2,r2  ’ 


while  satisfies 


(20a)  <P^^^'(s) 


(20b) 


1.  1.1 


(i)p(i)  ^  (i)g(i)  (i) 

^2  ^2,  ^2  ^2  1  ^  * 


(i),  (i-1).  ,  (i) 

<0  (3  )  =  U  . 


Equations  (19a)  and  (20a)  constitute  a  (coupled)  system  of  matrix  Ricatti 
equations. 

Of  course,  it  may  happen  that  the  solution  to  (19),  (20)  does  not 
exist  on  [s^^'^^.Sp].  It  will,  however,  exist  on  some  maximal  interval 
Cs^^^'^^t),  t  >  We  do  not  expect  to  continue  the  numerical 

integration  of  (19),  (20)  all  the  way  to  t.  Indeed,  let  A  >  1  be 


an  a  priori  given  constant  and  set  d^^^(3)  »  |<t>  '(s)j 


(p,n) 


|ld  ^(3)| 


(p,n) 


Let  t^^^  =  supls:  d^^^(s)  S  Ad^^^s^^  ^^)}  and 


then  set  s^^^  »  min{s2.t^^^}.  We  propagate  the  solution  of  (19),  (20) 
only  to  3  =  s^^^,  at  which  point  we  set 


e.  ^ 

.\v  .'-’N 
.>  .y. 


This  completes  Step  i. 


It  is  in  this  manner  that  we  traverse  the  interval  [31,32]  from 

left  to  right,  adaptively  constructing  the  (Riccati)  factorization  in  such 

a  way  that  it  remains  bounded  in  the  sense  of  Definition  2.14-2.  Indeed, 

T  ”1 

we  easily  have  that  in  the  spectral  norm  |(  $(s)  (3 ) )  |  ^  1  so  that 

|(<t(3)<fr'^(3))  ,  S  — 7— r  while  |'f^(s)|.  .  $  M  =  Ad,  where 

'  '(p.p)  ®(p)  '  '(:j),n) 


d  =  3up{|[Id  !  <»j| 


-  1  ‘fjlfn  :  I*  is  of  size  (p,n-p)  and  S  1). 

Pl'Pf”/  2-J 
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V  ^ 


We  refer  to  the  points 


,  i  «  1 , 2 , . . .  as  switching  or  resetting 

points. 

Remarks . 

1.  The  factorization  (15),  (16)  with  the  conditioning  matrix 
given  in  (17)  again  comprise  a  system  of  matrix  ODE's  on  a  manifold.  How¬ 
ever,  the  manifold  is  given  explicitly  by  (18).  The  advantages  from  a 
computational  point  of  view  are  striking: 

2 

a.  We  are  spared  the  integration  of  n^  equations  in  the 

2 

forward  direction  and  equations  in  the  backward 

direction. 

b.  The  constraint  manifold  (18)  is  maintained  explicitly. 

2.  Let  us  emphasize  the  very  important  point  that  the  perturbation 

A^(s)  of  Section  2.4  is  of  the  type  [0^  I  A].  That  is,  we  do  not  have 
perturbations  in  the  first  part  of  because  we  explicitly  maintain 

in  the  form  Cldp,<t].  This  point  is  not  exploited  by  the  general  analysis 
which  led  to  Theorem  4.1. 

3.  In  connection  with  the  Riccati  equation  approach  we  would  like 
to  mention  the  important,  but  generally  unknown,  papers  by  J.  Taufer.  See 
[22],  [23]. 

We  must,  on  the  other  hand,  confront  the  difficulties  associated 
with  the  numerical  integration  of  matrix  Riccati  equations.  We  have 
already  showed  how  to  exploit  the  possibility  that  the  Riccati  solution 
trajectories  may  "blow  up"  in  finite  time  (i.e.,  by  the  use  of  our  switch¬ 
ing  strategy).  We  turn  in  the  next  section  to  the  development  of  a  stable 
implicit  solver  for  Initial  value  problems  for  matrix  Ricatti  equations 
which  exploits  the  special  structure  of  the  quadratic  right  hand  side. 
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COMPUTATIONAL  ASPECTS  OF  THE  FACTORIZATION  METHOD 


Any  numerical  realization  of  the  method  of  factorization  should  be 
robust,  efficient,  accurate,  and  stable.  Moreover,  it  should  exploit  and 
preserve  the  special  character  of  the  particular  factorization  being  used. 
These  requirements  pose  special  problems  for  the  designer  of  factorization- 
based  codes. 

We  have  already  noted  the  advantages  of  explicit  over  implicit 
constraint  manifolds.  There  is  also  the  question  of  the  matrix  character 
of  the  initial  value  problems  of  the  factorization.  It  is  tempting  to  use 
one  of  the  excellent  modern  adaptive  initial  value  codes  now  available. 
This,  however,  means  that  the  matrices  must  be  "unrolled”  into  equivalent 
vector  form.  This  poses  no  special  problems  except  in  cases 
— usually  the  ones  of  practical  interest — in  which  a  stiff  (implicit) 
solver  is  required.  Propagation  of  size  (njjj,n)  transition  matrix  in 
unrolled  form  with  a  stiff  solver  requires  the  computation  and  repeated 
decomposition  of  Jacobian  matrices  of  size  (n^n,n^n).  This  represent  a 
severe  computational  burden  even  for  moderate  n^^  and  n  (n^^  =»  10  and 
n  =  20,  say) . 

Indeed,  this  computational  burden  is  so  severe  that  it  is  part  of 
the  folklore  of  control  engineering — where  matrix  Riccati  equations  play  a 
central  role — that  one  should  not  integrate  the  matrix  Riccati  initial 
value  problem  in  order  to  determine  steady-state  solutions.  Special 
methods  have  been  developed  to  solve  the  algebraic  (steady-state)  Riccati 
equation  instead  ([14],  [21]). 

We  will  show  in  this  section  how  these  computational  limitations 
can  be  overcome  by  the  design  of  special  matrix  initial  value  solvers 
which  exploit  the  structure  of  the  factorization  equations.  We  analyze, 


in  particular,  the  accuracy  and  numerical  stability  of  a  solver  for  matrix 


initial  value  problems  of  Riccati  type.  We  will  need  the  concept  of  a 
numerical  process  (cf.  [4]). 


3.1  Ntmerical  processes 

Denote  by  jjij  the  vector  space  of  matrices  of  size  (n,m) 

endowed  with  a  norm  I*),  and  consider  the  matrix  initial  value  problem 


(la)  4>'(3)  -  R(4>(3),3), 

(lb)  $(  Sq  )  =  iI'q  , 


with  [a.b]  5  3  -  <t(s)  €  and  R:  x  [a.b]  -  a 

Lipschitz  continuous  mapping. 

Definition  1.  A  one  step  numerical  process  in  is  a  sequence 

c  and  a  sequence  {hj^}  c  R  such  that 

(2) 

where  <tQ  and  the  sequence  {hj^}  are  given,  and  where  each  is  a 

mapping 


Pi'  ml  *  H  B 

1  (n,mj  ^n,m;  " 

We  regard  the  process  (2)  as  a  discrete  approximation  to  the  con¬ 
tinuous  process(l)  if  it  is  consistent;  that  is. 

Definition  2.  The  numerical  process  (2)  is  consistent  to  order  p  with 


(3a) 


*  h^)  +  0(h^  ’) 


^0  ^ 

J=0  J 


and  <ti(s)  satisfies  (la)  with  initial  condition  it(s^)  =  If  in  addi¬ 


tion  the  form  (2)  is  regular,  i.e., 


|Pi(<ti.hi)  -  P^(ii.hi)l  ^ 

for  ♦.  €  M„  _  with  |<t.  -  i,  |  <  e  for  all  i  and  with  M  independent 

1  ii  f  III  *1  X  ■ 

of  i,  then  (2)  is  a  general  one  step  method  for  the  solution  of  (1a,b) 

which  converges  with  the  rate  0(hP)  (cf.  [^],  Section  3.3.1).  B 

Another  way  to  say  this  is  that  (1)  is  the  (order  p)  closure  [5] 
of  the  numerical  process  (2). 

Definition  3.  The  one  step  numerical  process  (2)  is  bounded  if  there 


exists  M  ^  0  such  that 


^  M 


and  finally, 


Definition  4.  The  one  step  numerical  process  (2)  is  stable  if  there 
exists  L  i  1  and  eg  >  0  such  that 

(5)  i  L|$.  - 

whenever 


*1+1  “ 


If  L  <  1  then  (2)  is  strongly  stable. 


Definition  4  is  an  adaptation  of  the  notion  of  BN-stability  [9].  In 
contrast  to  the  concepts  of  A-stability  and  B-stability,  the  definition 
is  made  without  reference  to  a  particular  test  equation. 


3.2  A  one-step  process  for  Ricattl  equations 

Consider  the  special  case  of  (1)  given  by  (of.  equations  (2.5-29)) 


(la)  yUs)  -  -D(3)  +  A(3)y(s)  -  'i'(3)B(3)  +  'i'(3)C(3)4'(3) 


(1b)  'f(3Q)  -  fQ. 


Here  'l'(s)  €  A(s)  €  B(s)  €  C(3)  6 

and  D(s)  €  Let  us  use  the  notation  «  A(3j^  +  cxhj^),  for 

a  €  [0,1].  Then  we  define  a  one  step  numerical  process  for  (1)  in  two 


stages : 


Cld_  -  V,  h,(A.  „  +  ,  )] 


i'  i+  ‘4  i  i+  \ 


■f’l-  .4*  Vy ’’ 


-Cyi.  4  •  y.  »i.  y.-Oi.  y,’’ 

y  -  Cl.  y  ’1.  y  >]■'• 


Remarks. 


1.  It  is  certainly  possible  to  eliminate  the  intermediate  solu¬ 


tion  4' 


i+ 


between  (2a)  and  (2b)  in  order  to  write  (2)  in  the  form 


(3.1-2). 


2.  We  evaluate  the  matrix  coefficients  A,  B,  C,  D  only  once,  at 

the  midpoint,  s  .  of  the  step.  The  RHS  of  (1)  is  not  evaluated  at 
i+  \ 

all. 

3.  There  are  two  matrix  decompositions  per  step;  one  of  type 
(n,n)  and  the  other  of  type  (m,m),  for  an  operation  count  of 

0(n^  +  m^).  This  is  very  favorable  compared  to  the  0(nm)^  operations 
required  for  the  factorization  of  the  Jacobian  of  the  RHS  of  (la)  when 
written  in  unrolled  (vector)  form. 

Of  course,  in  an  implementation  of  the  Riccati  factorization  in 
order  to  solve  BVP's  we  must  propagate  a  transition  vector,  i</,  as  well 
as  the  transition  matrix,  'V.  In  addition  to  (1),  we  have  (cf.  equations 


2.5-20) 

(3a) 

lip'(s) 

-F(s)  A(s)i)»(3) 

-  't'(s)G(s)  +  4'(s)C(s)\t)(s) , 

(3b) 

>1^(30)  = 

'^0* 

in  which 

4^(3) 

,  F(s)  €  j  ) , 

G{3)  €  and  A  and  C 

are  as  above.  The  integration  of  (3)  is  done  simultaneously  with  that  of 
( 1 )  by  observing  that  (1),  (3)  can  be  written 


[4'  I  >1^] 


-CD  I  F]  +  hL'f  I  i|i] 

I  ' 


B  I  G  C 

-Ct  !  '!<]  *  ['?  j  - let  1  1)^] 

'  0  '  0  '  0  ' 


which  is  again  of  the  type  (1).  We  apply  the  scheme  (2)  to  the  augmented 


system  (U)  in  practice.  Of  course,  the  zero  coefficients  are  not  stored 
in  the  computer  implementation  of  (2)  applied  to  (i<). 


Remark. 

^  ^  ”(n,p)  °  ^  ”(m,p)'  ^  ”(n,p)  equations 

(3).  This  corresponds  to  solving  the  original  BVP  (2.1-1)  for  p  right 
hand  sides.  This  can  be  accomodated  trivially  in  (2)  (4). 


3«3  Consistency  of  the  method 

The  numerical  process  (3-2-2)  is  consistent  to  order  2  for 
(3.2-1).  Indeed  (3.2-1)  is  of  the  type 


(1) 


f(s; '?,'?), 


-^(30) 


while  the  corresponding  discrete  process  (3-2-2)  is  of  the  type 


(2a) 


i*  \ 


\  ^  .4  ‘’r’i.  4  > 


(2b) 


'i+1 


Vv,  y  -Vy,’' 


In  order  to  analyze  the  one-step  error  of  (2),  it  is  sufficient  to 
consider  the  case  i  »  0;  and  for  rotational  clarity  we  write  h  «  hQ 
and  suppress  the  dependence  upon  3  in  (1)  and  (2).  The  exact 

‘4 

solution  satisfies 


(3) 


4'(sQ+h)  -  'J'q  +  h*?^  +  ^  'J'q  O(h^) 


where 


WV"  «C"  TT 


I 


»• 

v'* 

i 

f- 

K- 


h.  • 

Pt 

.  • 
•/ 

t 


(>ia) 

^6 

(4b) 

■  e,<’o-''o>'’'’'o'''o' 

*  82<''o’''o'f<’o’’'o' 

(4c) 

o 

o 

■ 

•& 

(4d) 

g2 

On  the 

other  hand, 

from  (2)  we  have  that 

(5) 

= 

^0  *  1  >4  ^  " 

f('l', ,’i'  ,,  )] 

4 

Thus 


(6)  fCsQ+h)  -  H-, 


However, 


f(»  .» ) 


’o  ‘  I  *|82<’'o’’'o>'’'*o'’y,  > 

+  O(h^). 


«(»„,»„)  -IfOo.^t,)  -IfCf  .»,  ) 


•  r  i:8i<’o-''o>'''’o’’o>  *  e2<’'o'’o"'"o’’'o” 


r  U,(1  ,  .1  „  y  )  •  82(»o.»o)r(»o’''  y  ” 


+  O(h^), 


U3 


■ 'vS^ 

'■.*  ‘.-y 


+  O(h^) 


Therefore,  we  have 


(7)  '('(sQ-h)  -  ^  [g^('¥Q,'i'Q)f(4.o.'¥  ^  g2(’yo.4'Q)f('l'Q.'l'  )] 


r  V  y  ^^^'*'1’’*'  y  ^  ^  ^2%'^  y 


+  O(h^). 


But  also,  we  note  that 


,4'  ) 


f(4'Q,'l'Q)  +  0(h), 

g^('FQ,4'Q}  +  0(h), 

f(¥Q,4'Q)  +  0(h), 


SO  that  up  to  terms  of  order  0(h),  the  bracketed  quantities  in  (7)  sum 
to  zero; 


(8) 


H'(sQ+h)  - 


O(h^), 


In  practice,  the  2—  order  process  (3.2-2)  can  be  extrapolated  to 
higher  order,  with  the  order  and  the  step  size  h  varied  adaptively  based 
upon  comparisons  of  consecutive  stages  in  the  extrapolation.  An 
implementation  at  the  University  of  Maryland  uses  a  fixed  number  of 
extrapolation  stages,  three;  it  is  a  method  of  type  with  adaptive  step 
selection. 

The  step  control  is  based  on  the  usual  simultaneous  usage  of 


methods  of  two  different  orders.  The  local  error  Is  controlled  so  that 
the  method  uses  the  minimal  number  of  steps  when  the  perturbation  is 
measured  in  the  norm  L,  or  on  the  norm  (see  [7]).  In  [7]  it  has 

also  been  shown  that  the  error-per-unit  step  approach  is  optimal  with 
respect  to  the  L*  norm  and  the  error  per-step  is  optimal  with  respect  to 
the  norm.  We  give  some  numerical  examples  in  Section  U.  But  first 

let  us  analyze  the  stability  of  the  numerical  processes  (3. 2-2, 4). 

3>4  Stability  of  the  method 

The  boundedness  and  stability  of  the  process  (3.2-2)  are  inherited 
from  the  continuous  processes.  If  the  solution  to  (3.2-1)  exists  on 
Cs^,S2]  5  3o»  then  it  is  uniformly  bounded  on  [s^,S2].  Suppose  further 
that  the  solution  is  locally  stable  forward  in  s.  Then  it  is  not  hard  to 
see  that 

Lemma  1.  If  s  -»  y(s)  satisfies  (1)  on  Csi,S2]  5  Sq,  then  the  trajec¬ 
tory  'Vis)  is  locally  stable  forward  in  s  if  and  only  if 

(1 )  Re  X(a)  -  Re  w(s)  ^  0 

where  Xis)  is  siny  eigenvalue  of  'F(s)  *  A(3)  +  ’{'(3)0(3)  and  ii(3)  is 
any  eigenvalue  of  G(s)  =  B(s)  -  C(s)'?(s). 

Proof.  We  need  only  observe  that  (1)  is  necessary  and  sufficient  for  the 
trivial  solution  of  the  linearized  problem, 

(2)  <f'(3)  -  [A(s)  +  '1'(3)C(3)]<K3)  -  <K3)[B(3)  -  C(3)'l'(3)], 

to  be  locally  stable  forward  in  3. 

Indeed,  writing  (2)  in  unrolled  vector  form  we  analyze  the  eigen- 


values  of  the  Kronecker  sum  of  the  matrices  F(3)  and  G(s)  (see  [8],  p. 
230).  The  eigenvalues  of  this  sum  take  exactly  the  form  (1).  H 

Corollary  2.  The  trajectory  is  locally  stable  backward  in  s  if  and 
only  if 


Re  1(3)  -  Re  4(3)  i  0. 


This  stability  property  is  inherited  by  the  discrete  process 


(3.2-2).  Indeed,  let 


A  -  f 

*i  i  i 


where  solves  (3.2-2)  with  initial  value  Then  we  have 


a  A  -f  — fA  it  -tB  -t-tC  f 

‘4  i  2  ‘4  if  ‘4  i  if  ‘4  i  if  \  if  \ 


*1.  y.  S-  ^  ^  ’ 


ri . 

K  »A  f— i-fA  t  -tB  C  t. 

if  \  2  ^^f  y,  if  ‘4  ^fi“if  \  if  ‘4  ''if  \  ’if  4 


'^i-'-l^if  \  ^if  \ 


Upon  eliminating  we  have 


r  <‘i.  .4  *’i.  .4  =1.  4  -  r<‘i.  4  -’i.  4  ‘=1.  !4  • 


[Id  -  (3  -  C  ,  f  )][Id  +  ^  (B  ,  -  C  ,,4'^  .  )] 

1  m2  if  ‘4  if  ’4  4  ^  I*  4  4  4 

Now  set 


-  -  w  *  .  •  -%  A  -  •  *  ■  -  •  ^  •  h>  -  «  .T-  1  ^  - 


y  IT  II  w.y  Ml  ij  ijfiji  lum  »j  'j  pj  w  fj  pj  ^ 


*  A-  1/  *  '•'•  I,  C.  ,,  , 

i  i+  %  i+  %  i+  V, 


and  rewrite  (6)  as 


Vj  "  ^i+  V*  ’*'i+  \ 


h.  h  -1  h  h.  -1 

(8)  $  =  (Id  +^F)(Id  -^F)  $  (Id  -  ^  G,)(Id  +  ^  G  ) 

1  +  1  n^ii  nzi  iDiisl  ni£i 


••>y 

.'3 

>.•  >Vj 

•* V  A  ^ 

'  ‘Wl 


We  therefore  have  that 


(9) 


where 


hi., I  s  rJ*J, 


f  .*  .* 

c::!^ 


(10a) 


»  maxjYj, 


and 


(10b) 


i 

(1  +  ^  \^)  (1 


h. 

r 


h,  h . 

(1  -  (1  ♦  U,) 


In  (10),  the  maximum  is  taken  over  of  the  form  (106)  where  is 

any  eigenvalue  of  Fj^  and  is  any  eigenvalue  of  Gj^. 

For  numerical  stablity  of  (3-2-2)  we  must  ensure  that  |y^|  SI 

2 

(|Yj^|  <  1  for  strong  stability).  Let  us  compute,  then,  (Y^|  .  We  drop 
the  index  i  for  rotational  clarity. 


(11)  |y|' 


(1  +  I  A)(1  -  I  A)  (1  -  I  u)(1  -  I  ^) 
(1  -  I  A)(1  -  I  Y)  (1  +1  u)(1  +  I  u) 


h  ,h, 

2  P^2^ 

^  h  .h. 
2  ^^^2^ 


Here 


’J"  ' 


'**•'*• 

-  H  ■  3! 
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p(h)  -  2  Re(A-u)  +  h((X|  +  |p|“  -  4  Re  X  Re  y) 


+  2h^(|u|^Re  X  -  |X)^  Re  u)  +  h^lxl^lu)^ 


q(h)  -  2  Re(w-x)  +  h(lx|^  +  \v\^  -  4  Re  X  Re  u) 


+  2  h^(|x|^  Re  u  -  liil^  Re  X)  +  h^lxj^UI^. 


Stability  is  guaranteed  whenever 


h  -hv  h  r  /hv  /h.,  .  ^ 

2  ''  2^  *  2  ^  0 


with  strict  inequality  for  strong  stability.  Of  course,  the  represen¬ 
tation 


r(|)  *  4  Re(X-y)  +  h^(iuj^  Re  X  -  |x|^  Re  u) 


gives  us 


Leoima  3-  If  and  are  discrete  trajectories  satisfying  (3.2-2) 

and  such  that  Re  X^  ^  0,  and  Re  ^  0  where  Xj^  is  an  eigenvalue  of 
A  I  C  ,,  and  u.  is  an  eigenvalue  of  B,.  1,  -  C  ?  ,  , 

1+  ‘4  i+  ‘4  1+  ‘4  i  4  i+  i+  \ 

then  the  process  (3-2-2)  is  stable  forward  (hj^  >0).  If  Re  Xj^  S  0  and 
Re  Uj  ^  0  then  the  process  (3.2-2)  is  stable  backward  (h^^  <  0). 

Lemma  4.  Let  X^  and  be  as  in  Lemma  1 ,  but  with  the  weaker  hypothe¬ 

sis  Re(X.-u<)  ^  0.  Then  the  process  (3-2-2)  is  stable  forward  (h,  >  0) 


(16) 


Remarks . 

1.  Lemma  3  is  the  matrix  ODE  analog  of  A-stabillty  in  the 

constant  coefficient  linear  case  (CCs)  =  0). 

2  2 

2.  If  |u|  Re  A  -  I A I  Re  u  >  0  while  Re  A  -  Re  p  <  0  in  Lemma 
4,  it  is  interesting  to  examine  how  h^g^  depends  on  A.  For  example, 

2  -  y 

suppose  A  =»  ap  with  0  <  a  <  1.  Then  o  ^  .  Thus  as  long 

as  the  stronger  hypotheses  of  Lemma  1  are  not  violated  too  severely  (a 
near  1)  the  stability  limit  on  the  step-size  will  still  be  quite 

large. 

3.  For  a  well-posed  elliptic  boundary  value  problem  it  can  be  shown 

that  Re(e.v.  F)  ^  0  and  Re(e.v.  G)  i  0  for  the  forward  process  and 

Re(e.v.  F)  ^  0  and  Re(e.v.  G)  $  0  for  the  backward  process  .  .  .  exactly 

the  conditions  of  Lemma  1.  This  will  also  be  the  case  for  the  Riccati 

equations  which  arise  in  classical  linear-quadratic  optimal  control.  These 

T 

are  usually  posed  backward  and  have  the  additional  property  that  -G 
=  F  >  0. 
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4 


NUMERICAL  EXAMPLES 


In  this  section  we  give  several  examples  of  the  performance  of  a 
factorization  based  two  point  boundary  value  code  based  on  the  ideas 
presented  above.  We  have  chosen  some  examples  to  illustrate  the  robustness 
and  effectiveness  of  the  solver  on  problems  stemming  from  engineering.  We 
include  also  an  unstable  turning  point  problem  for  which  our  method  fails. 

We  have  already  noted  how  the  independent  and  parallel  structure  of 
the  forward  and  backward  factorizations  coupled  with  the  solution  of  local 
linear  systems  at  the  target  points  combine  to  minimize  storage  and  the 
computational  burden  of  the  method.  The  adaptive  mesh  selection  is  based 
on  a  single  solution  of  the  problem,  in  contrast  to  the  multiple-pass 
approach  to  mesh  refinement  used  in  global  methods  such  as  finite  differ¬ 
ences,  collocation,  or  finite  elements.  This  feature  also  greatly  reduces 
the  computational  costs. 

The  numerical  examples  illustrate  these  features,  but  highlight  the 
performance  of  the  method  on  problems  having  a  singular  perturbation  char¬ 
acter,  i.e.,  problems  the  solutions  of  which  exhibit  boundary  or  interior 
layers.  We  show  that  such  problems  can  be  solved  effectively  without 
special  handling  such  as  upwlnding  or  asymptotic  expansions. 

4.1  A  stable  singular  perturbation  problem 

Consider  the  problem 


(la) 

eu"(3)  +  u'(s) 

(1b) 

u(0)  -  ud) 

the  solution  of  which  is 


1-exp(  f-) 

(2)  u(3)  =  3  -  - -. —  . 

l-exp(-  ^) 

There  13  an  0(e)  boundary  layer  at  3=0. 

There  are  a  number  of  ways  that  ( 1 )  can  be  cast  into  the  first 
order  form  2.1-1;  we  explore  some  of  the  theoretical  and  computational 
implications  of  such  re-formulations. 

Method  1  (Ml) 

T 

Let  w  =  (w^,W2)  be  defined  by 

(3)  w^  =  u,  W2  =  u' 
and  obtain 


"o  1  ■ 

"0  ■ 

(^a,b) 

B  » 

0  -  i 

• 

F  = 

_  1_ 

L  e  J 

.  e  _ 

(Uc.d) 

Ui  =  U2 

=  Cl 

0], 

Ui  = 

U2 

Method  2  (M2) 

T 

Let  w  =  (w^,W2)  be  defined  by 


(5) 

W-|  =  u, 

W2  -  EU ' 

and  obtain 

'o  i' 

e 

0 

(6a, b) 

B  = 

0  -i 

e 

,  F  - 

-1 

"  m 

- 

(6c,d) 

Ui  =  U2 

»  [1 

0],  Ui  - 

U2 
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(Ml)  and  (M2)  have  the  form 


(7a) 

(7b) 


w 


'  Bw  -  F, 
U-,w(0)  -  U2W(1)  -  0. 


The  two  possibilities  for  the  Riccati  factorization  (which  are  adaptively 
alternated)  corresponding  to  (7)  are  factorization  FI  : 


(8a) 

*1  = 

Id, 

(8b) 

A  ^  s 

2 

-B  +  B 

12  “lV2 

*2®22  *2®21*2' 

(8c) 

(p'  =. 

-^1  *  - 

*^2^2  *2®21  ’ 

and  factorization  F2: 

(9a) 

“®21  "  ®22^ 

■  *1^1  *  *1®12*1' 

(9b) 

$2  - 

Id, 

(9c) 

<p'  “ 

-^2*^22  - 

*1^1  ^  *1^12  ' 

where 

$  -  C'l'^  *2^ 

<p  are  the 

transition  matrix  and 

vector 

,  respectively. 

Suppose  that  matrix  B  and  vector  F  are  perturbed  by  b  and 
f,  respectively.  Then  the  solution  w  to  (7)  is  perturbed  by  v 
satisfying 


(10a) 

(10b) 


Utv(O) 


0. 


V* 


Bv  +  (bw  -  f) 


The  matrix  b  and  the  vector  f  arise  computationally  due  to  the 
discretization  errors  involved  in  solving  (8)  or  (9).  In  the  case  of  FI, 


if  the  discretization  error  of  (8b)  is  ~i>i2  of  (8c)  is  -f ■( , 

then  we  see  that 


(11a) 


(11b) 


For  the  factorization  F2  the  corresponding  perturbations  are 


(12a) 


0 

0 


(12b) 


where  now  -b2i  is  the  discretization  error  of  (9a)  and  -f2  is  the 
discretization  error  of  (9c).  (Eqns  (8a)  and  (9b)  are  ’’solved  "  exactly.) 

Since  for  Ml,  '^2  ‘  boundary  layer,  we  see  that  bi2 

=>  0(A)  perturbations  in  8^2  lead  to  large  0(A/e)  errors  in  v.  We 

therefore  would  like  to  ensure  that  9^2  ”  0(Ae)  in  the  layer  in  order  to 
get  V  »  0(A).  This  can  be  done  by  brute  force  using  an  integration  toler¬ 
ance  T  »  Ae.  This  strategy  will  cause  the  adaptive  solver  to  use  more 


-\-v' 

y.‘-v 


V 
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steps  globally  as  well  as  In  the  boundary  layer.  The  problem  does  not 
arise  in  factorization  F2  since  the  perturbation  has  been  shifted  to  b2i i 
the  coefficient  of  -  0(1)  in  (10a).  Moreover,  although  the  factori¬ 

zation  FI  is  in  effect  at  the  start  of  both  forward  amd  backward  sweeps  due 
to  the  structure  of  the  boundary  condition  matrices  and  U2,  the 

rapid  growth  of  *2  (for  e  <<  1 )  causes  a  switch  to  F2  after  only  a  few 
steps. 

Now  consider  the  effect  of  measuring  the  error  by  the  norms 

(13a)  lv|^  -  sup{[v^(s)l  +  1v2(3)1}, 

s 

(13b)  |v|2  -  sup{|v^(s)|  +  e|v2(3)|}. 

s 

The3e  vector  norms  induce  corresponding  matrix  norms  on  b  given  by 
|b|^  -  sup  maxC|b^  J  Cb^2l 

3 

(14b)  |b|2  »  3upmaxC|b,J  +e|b2j,  e'^|bi2l  l^22l^’ 

s 

For  factorization  FI,  keeping  |b|2  <  A  means  exactly  that 
|b^2l  ^  Consider  a  transformed  problem 

(15)  w'  -  Bw  -  F 

for  which  |b|^  »  I^l2-  ^  ’  ^0  e^’  ^  ^  * 

ABA“^ ,and  F  «  AF.  Then  it  is  easy  to  see  that  if  (7)  is  solved  with  B,  F 

replacing  B  and  F,  then  b^2  *  ^  ^^12*  *  |b|i  - 

sup|b^2l  *  3up|e  solving  FI  with  tolerances  A,  5 

S  3 

for  (15)  is  equivalent  to  solving  (7)  with  tolerances  eA,  e6.  Of  course 


B  and  F  lead  exactly  to  problem  M2. 

We  have  shown,  then,  that  transformations  of  the  original  problem 

(1)  into  various  I®''  order  forms  can  be  interpreted  as  selecting  a  norm  for 

the  original  problem,  and  affect  the  interpretation  of  the  perturbations  of 

the  coefficients  caused  by  discretization  errors  in  the  solution  of  the 

factorization  initial  value  problems. 

We  further  consider  four  other  first  order  problems  computing  w  - 
T 

(w^ ,^2)  : 

Method  3  (M3) 

(16)  =>  eu'  +  u,  W2  =  u' ; 

Method  4  (M4) 

(17)  =■  u,  W2  =•  cu'  “  u; 

Method  5  (M5) 

(18)  “  u,  '^2  “ 

Method  6  (M6) 

(19)  »  u,  '^2  °  +  (3  -  ‘4  )• 

We  are  interested  in  the  computation  of  u'(0)  by  the  mentioned  six 
formulations.  Therefore  only  one  target  point  is  considered,  namely  s  = 

0.  Hence,  only  integration  from  right  to  left  has  to  be  performed  by  the 
factorization  method.  The  computation  has  been  in  double  precision  and 
"per  step"  step  selection  criterion.  The  initial  step  was  taken  to  be  h 
»  e.  Table  4. 1-4.6  show  the  error  in  u’(0)  obtained  by  the  method  (M1)- 
(M-6)  for  various  e  and  tolerances  t.  In  addition,  the  number  of 


integrations  steps  is  given  in  the  tables. 

Comparing  Tables  4.1  and  4.2  we  see  that  the  methods  (Ml)  and  (M2) 
produce  identical  results  if  -  £^2  where  by  the  tolerances  used 

in  method  Mj  were  denoted.  This  is  directly  related  to  the  analysis 
mentioned  above.  Methods  (Ml)  and  (M3)  are  using  as  one  of  the 
variables  W2  -  u',  while  (M2),  (M4),  (M5),  (m6)  are  using  the  variable 
W2  involving  eu'.  This  leads  to  similar  performances  of  these  two 
groups  of  methods.  Nevertheless,  there  are  some  differences  (see  (M2)  and 
(m6))  which  are  caused  by  the  different  structure  of  the  perturbations  in 
B  and  F.  Results  shown  in  Table  4.1  are  heavily  influenced  by  low 
tolerances  x  for  which  the  number  of  steps  is  independent  of  t.  This 
disappears  for  t  smaller  as  can  be  seen  in  Table  4.2,  which  is,  as  we 
said,  essentially  the  method  (Ml)  with  tolerance  xe.  We  see  in  Table  4.2 
that  the  error  in  u'(0)  for  fixed  c  is  proportional  to  x.  This  is 
because  the  perturbations  in  the  input  data  caused  by  the  approximate 
solution  of  the  solved  ODEs  are  of  magnitude  x. 

From  (10)  we  thus  expect  that  the  error  in  u'(0)  is  of  order 
We  see  this  character  from  Table  4.2  especially  for  small  x  and  e 
(when  we  are  in  the  asymptotic  range).  For  a  particular  value  of  x  (x  • 
10  °)  there  is  a  (local)  increase  in  accuracy.  This  is  the  effect  of 
some  cancelation  and  was  observed  also  for  some  other  sequences  of 
tolerances.  The  magnitude  of  the  error  can  be  computed  from  (10)  assuming 

that  |b,,  I,  ,  If  I,  *  nx  where  n  is  the  number  of  steps.  This 

^1 

estimate  is  directly  related  to  the  fact  that  "per  step"  strategy  leads  to 
the  optimal  distribution  of  steps  minimizing  the  perturbation  in  the  Li~ 
norm.  See  [7]. 

The  methods  (M4),  (M5),  (M6)  show  similar  performance  although 
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method  (M5)  gives  better  results  than  (M4),  (M6).  The  factor  common  to 


(M4)  and  (M6)  is  the  relation  of  u  and  ^2  with  a  negative  sign  in  the 
differential  equation  (for  s  •  0  where  the  boundary  layer  is  located) 
while  (M5)  uses  the  relation  with  opposite  sign.  The  different  character 
of  perturbations  then  affects  the  error  in  u'(0).  We  see  that  the 
transformation  of  the  original  problem  into  the  first  system  and  the 
choice  of  norm  affects  the  error  in  u'(0)  because  it  leads  to  different 
characterizations  of  the  discretization  errors  as  perturbations  of  the 
input  data. 

Table  4.1 

Method  1.  Error  in  u'(0)  and  number  of  steps  backward. 

josolutt  tolerance 


epsilon 

i.ooooe-01 

l.OOOOE+00 

3.O02OE-O6 

1 

1.5190E-O4 

3 

l.OOOOE-01 

i.ooooe-02 

1.1043E-03 

4 

l.OOOOE-03 

l.a64lE-0S 

6 

i.ooooe-04 

a.3oiaE-oi 

7 

i.ooooe-05 

2.2905E+00 

9 

l.OOOOE-06 

2.2825E+01 

11 

i.ooooe-07 

2.2941E+02 

12 

l.OOOOE-02 

l.OOOOE-03 

3.0B20E-06 

3.0820E-06 

1 

1 

1.5190E-04 

1.7454e-05 

3 

4 

9.7234E-06 

5.4961E-06 

5 

5 

1.8641E-02 

3.4241E-03 

6 

E 

2.3018E-01 

2.30ia£-0l 

7 

7 

2.2905E+00 

2.2905E+00 

9 

9 

2.2KS+01 

2.2825E+01 

11 

11 

2,294l£+02 

2.294a£+02 

12 

12 

l.OOOOE-04 

l.OOOOE-05 

1.024flE-07 

3.94E8E-08 

2 

3 

2.0155E-05 

1.7975E-06 

6 

B 

2.129K-06 

6.0009E-08 

6 

8 

3.7770E-05 

2.0772E-05 

7 

8 

3.5308E-02 

1.9995E-03 

8 

8 

2.2905E+<» 

3.3971E-01 

9 

9 

2.2825E+01 

2.2S25£-^01 

11 

11 

2.2941E+08 

2.2941E+02 

12 

12 

l.OOOOE-06 

l.OOOOE-07 

l.2039e-O8 

1.7953E-09 

5 

7 

1.0358E-07 

i.aTsgE-oa 

13 

20 

3. 1869E-09 

2.0307E-09 

11 

17 

2.2579E-07 

9.4277E-0a 

10 

13 

7,3388E-04 

:.70oeE-o5 

9 

11 

2.0070E-02 

7.3869E-03 

10 

11 

3.3986£^ 

2.020eE-01 

11 

11 

2.2941E+02 

3.9009E+01 

12 

12 

7.4ii3E'02 

12 

2.0208Ef0O 


Table  4.2 

Method  2.  Error  in  u'(0)  and  number  of  steps  backward. 


l.OOOOE-01 

l.OOOOE-OE 

epsilon 

l.OOOOE+00 

3.0a20E-0& 

3.0e20E-06 

1 

1 

l.OOOOE-01 

1.5190E-04 

1.74S4E-0O 

3 

4 

i.ooooe-02 

S.4%lE-0fi 

2.1295E-06 

5 

a 

i.ooooe-03 

3.7770E-05 

a.077a£-05 

7 

a 

l.OOOOE-04 

1.999SE-03 

7.338aE-04 

a 

9 

l.OOOOC-05 

2.0070E-02 

7.3669E-03 

10 

11 

l.OOOOE-06 

3.020aE-01 

7.4113E-02 

11 

12 

i.ooooe-07 

2.0£OaE>00 

7.4116£-0l 

13 

14 

aasolutt  toliranct 


l.OOOOE-03 

l.OOOOE-04 

l.OOOOE-OS 

3.0a20E-06 

1.0a48E-07 

3.94aaE-08 

1 

2 

3 

2.0155E-05 

1.797S-0a 

1.03SaE-O7 

a 

8 

13 

a.0009E-08 

3.18£9E-09 

2.030aE-09 

a 

11 

17 

2.2579E-07 

9.4278E-08 

a.3E24E-08 

10 

13 

IB 

1.700aE-05 

i.iagaE-oa 

a.2e31E-07 

11 

14 

20 

1.7051E-04 

1.1849E-05 

a.4ooaE-oa 

13 

la 

21 

1.7221E-03 

1.1870E-04 

a.3983E-05 

14 

17 

23 

1.7222E-02 

l.iaa9E-03 

a.3994E-04 

la 

19 

24 

i.ooooE-oe 

l.OOOOE-07 

l.OOOOE-08 

1.2039E-0S 

1.7953E-09 

3.087aE-10 

S 

7 

11 

1.27S9E-08 

1.13S2E-10 

l.a3SEE-10 

20 

31 

SO 

1.8774E-10 

2.13iaE-13 

3.3111E-12 

25 

38 

59 

a.3a83E-09 

3.410eE-l2 

l.a64SE-ll 

2e 

39 

ai 

a.3024E-08 

3.ai99E-ll 

a.d394E-10 

2B 

41 

aa 

a.4434E-07 

3.0559E-10 

a.9B49E-09 

29 

43 

a4 

a.440iE-oa 

3.1432E-09 

7.0315E-08 

31 

44 

as 

a.443aE-05 

2.7940E-0a 

a.a54SE-07 

32 

4a 

a7 

Table  4.3 

Method  3.  Error  in  u'(0)  and  number  of  steps  backward.  , 


epsilon 

l.OOOOE-01 

l.OOOOC-02 

l.OOOOE+00 

1.1738E-05 

1.1738E-05 

1 

1 

I.OOOOE-01 

1.5190E-04 

1.5190e-04 

3 

3 

l.OOOOE-02 

1.1043E-03 

9.7234£-0a 

4 

5 

l.OOOOE-03 

i.aa4iE'02 

i.ae4iE-02 

a 

a 

t.OOOOE-04 

2.30iaE-01 

2.30ia£-0l 

7 

7 

l.OOOOE-05 

2.2905E+00 

2.2905E+00 

9 

9 

l.OOOOE-06 

2.2a2SE-H}l 

2.2a25E+01 

11 

11 

l.OOOOE-07 

2.2941E+02 

2.2941E+02 

12 

12 

absolute  toltrann 


l.OOOOE-03 

l.OOOOE-04 

l.OOOOE-OS 

i.07aeE-o7 

4.1489E-08 

8.a509E-09 

2 

4 

5 

1.7454E-05 

2.0155E-05 

1.7975E-06 

4 

6 

a 

s.49aiE-oa 

2.1295EH)a 

a.0009E-08 

5 

a 

a 

3.4241E-03 

3.7770e-05 

2.0772E-05 

a 

7 

a 

2.30iaE-01 

3.5308E-02 

1.9995E-03 

7 

a 

a 

2.2905E+00 

2.2905E+00 

3.a971E-01 

9 

9 

9 

2.2a25E+01 

2.2825E'KI1 

2.2a2SE-H)l 

11 

11 

11 

2.2941E+02 

2.2941E+02 

2.2941E-m)2 

12 

12 

12 

l.OOOOE-06 

l.OOOOE-07 

1.00002-08 

4.2592E-09 

a. 78392-10 

1.27822-10 

a 

9 

14 

1.0358E-07 

1.2759E-08 

1.13582-10 

13 

20 

31 

3.iaaaE-09 

2.0307E-09 

1.87842-10 

11 

17 

25 

2.2579E-07 

9.4278E-08 

8.38242-08 

10 

13 

18 

7.33882-04 

1.7008E-05 

1.10982-08 

9 

11 

14 

2.0070E-02 

7.38892-03 

1.70512-04 

10 

11 

13 

3.a98aE-f00 

2. 02082-01 

7.41132-02 

11 

11 

12 

2.2941E+02 

3.90092+01 

2.02082+00 

12 

12 

13 

Table  4.4 

Method  4.  Error  in  u'(0)  and  number  of  steps  backward. 

abioluta  tolnranci 


l.OOOOE-01 

i.ooooe-oe 

l.OOOOE-03 

l.OOOCe-04 

l.OOOOE-05 

l.OOOOE-06 

l.OOOOE-07 

epsilon 

1.0000E-K>0 

.2.10a2E-0£ 

2.10e2E-06 

2.10e2E-06 

S.1426E-0a 

4.ia99E-08 

1.3130E-08 

2.494a£-09 

1 

1 

1 

2 

3 

4 

7 

l.OOOOE-01 

1.4772E-04 

2.3294E-04 

4.9107E-0S 

8.4943E-06 

4.4914E-07 

S.3726E-0e 

4.6S13EH)9 

3 

3 

5 

7 

10 

16 

25 

l.OOOOE-02 

5.5036E-O5 

4.9424E-09 

7.6923E-07 

a.2960E-Oa 

1.1879E-09 

8.3691E-10 

3.0013E-11 

5 

5 

7 

9 

13 

20 

30 

l.OOOOE-03 

2.02S3E-03 

a.f70l£-05 

2.2672E-05 

1.789eE-06 

3.6016E-08 

2.8657E-09 

4.4145E-10 

6 

7 

a 

11 

15 

22 

32 

l.OOOOE-04 

2.0639e'O2 

1.1727E-03 

2.2109E-04 

2.7409E-05 

3.6452E-07 

1.9823E-07 

7, 1304E-20 

a 

a 

10 

12 

17 

23 

33 

I.OOOOE-05 

2.132SE-01 

1. 2697E-02 

2.2S62E-03 

2.7437E-04 

2.4974E-06 

1.6795E-06 

7.360a£-07 

9 

10 

11 

14 

la 

25 

35 

l.OOOOE-06 

2.1326£>00 

1.1732E-01 

2.2541E-02 

2.Sa64£-03 

1.2741E-04 

7.322aE-05 

1.5972E-05 

11 

11 

13 

15 

20 

26 

36 

l.OOOOe-07 

2.133SE^1 

1.1833£>00 

2.1307E-01 

2.6545£H)2 

2.aS26E-03 

a.6016E-4)4 

5.2899E-07 

12 

13 

15 

17 

21 

20 

38 

Table  4.5 

Method  5.  Error  in  u'(0)  and  number  of  steps 

1  backward 

• 

1.00006-01 

1.00006-02 

1.00006-03 

absolute  tolerance 
1.00006-04  1.00006-05 

1.00006-06 

1.00006h)7 

epsilon 

l.OOOOE+00 

1.17386-05 

1.17386-05 

1.07666-07 

4.1489E-08 

8.65096-09 

4.25926-09 

6.76396-10 

1 

1 

2 

4 

5 

6 

9 

l.OOOOE-01 

3.56466-05 

7. 55456-05 

2.81186-05 

2.33436-06 

1. 1942E-07 

9.62846-09 

9.82806-11 

5 

5 

7 

11 

15 

22 

34 

l.OOOOE-02 

1.1666E-0S 

3.0201E-06 

5.27906-08 

1.0821E-09 

4.36336-10 

3.85826-11 

6. 65076-12 

7 

7 

10 

14 

19 

27 

41 

l.OOOOE-03 

4.6703E-04 

4.01196-05 

1.54906-06 

4.7685E-0e 

1.23386-09 

9.94316-10 

1.71446-10 

a 

9 

11 

16 

21 

28 

42 

1.00006-04 

4.76956-03 

3.82666-04 

1.46056-05 

5.7743E-07 

1.07956-07 

9.35146-09 

1.62256-09 

10 

11 

13 

17 

22 

30 

44 

1. 00006-05 

4.9322E-02 

4.68686-03 

1.4541E-04 

5. 73676-06 

1.08176-06 

9.33366-06 

1.61966-08 

11 

12 

IS 

19 

24 

32 

46 

l.OOOOE-06 

4.9322E-01 

4.68836-02 

1.48206-03 

5.73606-05 

1.09356-05 

9.48326-07 

1.64386-07 

13 

14 

16 

21 

25 

33 

47 

l.OOOOE-07 

4.9335E+00 

4.6923E-01 

1.4820E-02 

5.7447E-04 

1.09356-04 

9.47906-06 

1.65406-06 

14 

25 

18 

22 

27 

35 

49 

T^PfodueeTTronT 


l.OOOOE-08 

4.n36£-20 

11 

1.7075E-10 

39 

3.3396E-12 

47 

6.ia46E-ll 

49 

1.5270E-0a 

50 

1.8963E-07 

ss 

2.1769E-05 

53 

2.947&E-03 

55 


l.OOOOE-08 

1.2762E-10 

14 

1.5913E-10 

54 

1.4211E-14 

64 

4.4793E-11 

65 

4.5657E-10 

67 

4.6712E-09 

66 

4.6915E-08 

70 

4.7311E-07 

71 
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Table  4.6 

Method  6.  Error  in  u'(0)  and  number  of  steps  backward. 

absolute  tolsranc* 


epsilon 

l.OOOOE-01 

1.00002-02 

1.00002-03 

1.00002-04 

1.00002-05 

1.00002-06 

1.00002-07 

1.00002-08 

l.OOOOE+00 

6.27%E-<S 

6.27922-05 

167942-06 

1.13082-08 

121302-08 

6.60432-09 

1.42412-09 

2.37922-10 

1 

1 

2 

3 

5 

7 

11 

17 

l.OOOOE-01 

2.047aE-<H 

2.06042-04 

2.21922-05 

2.14502-06 

1.79082-07 

5.23192-08 

2.25382-09 

4.22612-11 

3 

5 

6 

10 

15 

23 

36 

57 

l.OOOOE-Oe 

2.1tl6£-05 

1.27092-05 

2.70522-07 

6.18352-08 

1.34352-09 

6.46522-10 

7.43232-12 

9.80552-12 

6 

7 

9 

13 

19 

27 

42 

65 

l.OOOOE-03 

2.95a6E-04 

6.28532-05 

2.46162-05 

101312-07 

4.67172-08 

1.57892-08 

8.54932-11 

7.41242-11 

a 

a 

10 

14 

20 

28 

43 

65 

l.OOOOE-04 

3.7127E-03 

6.15802-04 

2.48632-04 

2.61672-06 

7.05092-07 

1.72812-07 

3.22322-09 

1.38432-09 

9 

10 

12 

16 

21 

30 

44 

67 

i.ooooe-<s 

3.69512-02 

5.66952-03 

2.55892-03 

150992-05 

7.17162-06 

2.01762-06 

1.26212-07 

2.08752-07 

11 

12 

13 

17 

23 

31 

46 

68 

l.OOOOE-06 

3.aa24£-oi 

6.2a07E-02 

2.55872-02 

145382-04 

7.84662-05 

1. 18182-06 

1.36312-06 

1.06752-05 

13 

13 

15 

19 

24 

33 

47 

70 

i.ooooe-07 

3.70112+00 

6.28792-01 

2.57542-01 

2.11432-03 

5.47212-04 

8.29712-04 

1.33912-03 

1.5071E-<'' 

14 

15 

16 

20 

26 

34 

49 

71 

4.2  An  unstable  singular  perturbation  problem 

Consider  the  singular  perturbation  problem  of  turning  point  type 

(la)  ew"  +  asw*  -  0,  -a  <  s  <  b, 

(lb)  _  w(-a)  -  1 ,  w(b)  -  2. 

For  a  «  1  the  problem  is  stable  and  has  an  interior  layer  (shock)  at  s 

»  0  (Fig.  4.1)  for  all  a,b  >  0.  The  solutions  shown  in  Fig.  4.1  were 

-4 

computed  with  the  tolerance  t  »  10  for  all  shown  e. 


^  r  • 


In  the  case  a  » -1  a  boundary  layer  occurs  at  s  -  -a  if  a  >  b  and  s 
-  b  if  a  <  b.  If  a  »  b,  then  there  is  a  boundary  layer  at  each  end. 
The  solution  in  this  case  (a  »  b)  is  unstable  with  respect  to  the  data 
as  e  -*•  0.  The  exact  solution  is  antisymmetric  with  respect  to  the 
value  1.5.  It  can  be  computed  for  s  >  0  by  solving  the  (stable) 
problem 

(2a)  ew"  -  sw'  -  0 

(2b)  w(0)  -  1.5,  w(1)  -  2. 

Fig.  11.2  shows  the  solution  for  various  e  computed  by  the  factorization 

_ji 

method  with  tolerance  t  »  10  .  Solving  the  original  problem  (1)  with 

a  =•  b  =  1 ,  one  has  to  expect  that  the  results  will  be  very  poor  for  e 
small  because  of  the  instability  of  the  problem  (1).  Fig.  4.3  shows  that 
in  fact  for  e  small  the  factorization  method  completely  fails.  This 
example  shows  that  the  stability  of  the  problem  is  a  necessary  condition 
for  the  factorization  method  to  give  high  quality  results.  This  condition 
is  directly  related  to  the  interpretation  of  the  numerical  solution  as  the 
exact  solution  of  a  perturbed  problem. 
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where  is  the  Legendre  polynomial  of  degree  k.  Symmetry  about  y  » 

0  in  the  problem  (1),  (2)  suggests  an  approximation  for  w  of  the  form 

N 

(5)  w(x,y)  =■  I  w  (x)J,-  (y) 

j=0  J 

for  (x,y)  n.  The  convergence  and  approximation  properties  of  this 
method,  called  dimensional  reduction,  are  discussed  by  Vogelius  and 
Babuska  [22].  The  relevance  of  (5)  for  our  purposes  is  that  it  reduces 
the  PDE  (1),(2)  to  a  system  of  ODE's  for  the  vector  function 

(6)  W(*)  3  (Wq(  •)  ,w^  (•),...  ,Wj^( 

mapping  [0,1]  into  .  The  function  W  is  the  solution  of  the  two 

point  boundary  value  problem 

2 

(7a)  -hA  h~^BW(x)  =  G(x),  x  €  (0,1) 

dx 

(7b)  W(0)  =  W(1)  =  0 

where  G(x)  =  ( 2g(x) , 0 , . . . , 0) and  the  (N+1)  x  (N+1)  matrices  A  and 
3  are  the  mass  and  stiffness  matrices,  respectively,  of  the  basis 
That  is. 


(8) 


We  note  that  the  presence  of  the  parameter  h  in  (7a)  gives  the  TPBVP  (7) 
a  singular  perturbation  character  if  h  <<  1 . 


Figures  ^.4  and  4.5  are  plots  of  the  solution  components 

{w,.}  and  {will  for  the  case  N  »  4  and  h  ■  0.5, 

k=0 . N  k-0 . N 

with  g  given  by 

11,  0.475  ^  X  S  0.525 

0,  otherwise. 


Only  the  solution  on  0  ^  x  ^  0.5  is  plotted  since  w(x,y)  in  this  case 
is  symmetric  about  x  -  0.5  and  y  -  0.  Evidently  the  component  Wg 
dominates  the  solution  except  in  the  interior  layer  caused  by  the  flux, 
g.  A  contour  plot  of  w  obtained  from  expansion  (5)  for  h  «  0.5  and 
N  =•  4  is  shown  in  Figure  4.6.  The  singularity  at  the  corner  due  to  the 
step  in  the  flux  is  evident. 

-4 

The  computations  were  done  with  tolerance  of  t  -  10  .  In  view 

of  (3.2-4),  this  amounts  to  computing  the  exact  solution  to  the  first 
order  system 


(11a) 


^  W^  (x)  -  [Idjj  +  b(x)]W2(x)  -  f(x) 


(llb)  -hA  ^  W,(x)  +  h'’BW,(x)  -  G(x) 

dx  2  1 

(llc)  Wi(0)  -  W^d)  -  0 

where  W^  ■  W  and  the  perturbing  matrix  b(x)  and  vector  f(x)  are 
unknown  but  satisfy 


(12) 


lb,.(x)  ,  |fi(x)|  <  T, 


i,J  -  0,N,  X  €  (0,1 ). 


This  in  turn  implies  that  the  numerical  solution  is  the  exact 

solution  for  problem  (1)  in  which  the  flux  G  is  replaced  by  the  flux 


G(x)  -  hA  ^  {bW^(x)  +  f(x)). 
dx  2 


unctions 
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